#####################################
# Install/load packages and functions
#####################################

# remove objects
rm(list=ls())
# detach all libraries
detachAllPackages <- function() {
  basic.packages <- c("package:stats", "package:graphics", "package:grDevices", "package:utils", "package:datasets", "package:methods", "package:base")
  package.list <- search()[ifelse(unlist(gregexpr("package:", search()))==1, TRUE, FALSE)]
  package.list <- setdiff(package.list, basic.packages)
  if (length(package.list)>0)  for (package in package.list) detach(package,  character.only=TRUE)
}
detachAllPackages()

# load libraries
pkgTest <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[,  "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg,  dependencies = TRUE)
  sapply(pkg,  require,  character.only = TRUE)
}

# here is where you load any necessary packages
# ex: stringr
# lapply(c("stringr"),  pkgTest)

lapply(c("car", "cregg", "data.table", "estimatr", "lmtest", "memisc", 
         "multcomp", "multiwayvcov", "pBrackets", "stringr", "texreg"),  pkgTest)

# set working directory
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# set seed in honor of T-Swift
set.seed(1989)

###################
# load and fix data
###################

data <- fread("funckmccabe_data.csv", header = TRUE, stringsAsFactors = FALSE)
data <- data[-c(1:2),]

# recoding respondent's placement on a seven-point PID scale as numeric
data$pid7 <- ifelse(data$Q13=="Strong Democrat", 7, NA)
data$pid7 <- ifelse(data$Q13=="Not a very strong Democrat", 6, data$pid7)
data$pid7 <- ifelse(data$Q15=="Democratic", 5, data$pid7)
data$pid7 <- ifelse(data$Q15=="Neither", 4, data$pid7)
data$pid7 <- ifelse(data$Q15=="Republican", 3, data$pid7)
data$pid7 <- ifelse(data$Q14=="Not a very strong Republican", 2, data$pid7)
data$pid7 <- ifelse(data$Q14=="Strong Republican", 1, data$pid7)

# coding whether respondents were assigned to the forced choice or abstention
# arms of the replication
data$exp_arm <- ifelse(data$FL_180_DO!="", "Forced", "Abstention")

# saving respondent-level information about party ID and experimental arm 
# for later use
demo_data <- subset(data, select = c(ResponseId, pid7, exp_arm))

# for each of the six tasks respondents completed, the following six code chunks
# use variables provided by Qualtrics to determine the partisanship of each
# profile presented, whether the task was situated in a high, moderate, or low
# information environment, and extracts the choice the respondent made in the task
data$party1_A <- ifelse(str_detect(data$FL_180_DO, "D-R") | str_detect(data$FL_186_DO, "D-R"), "Democrat", NA)
data$party1_A <- ifelse(str_detect(data$FL_180_DO, "R-D") | str_detect(data$FL_186_DO, "R-D"), "Republican", data$party1_A)
data$news1_A <- data$`G-1-1-1`
data$party1_B <- ifelse(str_detect(data$FL_180_DO, "D-R") | str_detect(data$FL_186_DO, "D-R"), "Republican", NA)
data$party1_B <- ifelse(str_detect(data$FL_180_DO, "R-D") | str_detect(data$FL_186_DO, "R-D"), "Democrat", data$party1_B)
data$news1_B <- data$`G-1-2-1`
data$info_env1 <- ifelse(str_detect(data$FL_180_DO, "High") | str_detect(data$FL_186_DO, "High"), "High", NA)
data$info_env1 <- ifelse(str_detect(data$FL_180_DO, "Moderate") | str_detect(data$FL_186_DO, "Moderate"), "Moderate", data$info_env1)
data$info_env1 <- ifelse(str_detect(data$FL_180_DO, "Low") | str_detect(data$FL_186_DO, "Low"), "Low", data$info_env1)
data$choice1 <- ifelse(data$exp_arm=="Forced", paste0(data$Q86, data$Q92, data$Q96, data$Q104, data$Q100, data$Q108),
                       paste0(data$Q235, data$Q239, data$Q243, data$Q247, data$Q251, data$Q255))

data$party2_A <- ifelse(str_detect(data$FL_181_DO, "D-R") | str_detect(data$FL_187_DO, "D-R"), "Democrat", NA)
data$party2_A <- ifelse(str_detect(data$FL_181_DO, "R-D") | str_detect(data$FL_187_DO, "R-D"), "Republican", data$party2_A)
data$news2_A <- data$`G-2-1-1`
data$party2_B <- ifelse(str_detect(data$FL_181_DO, "D-R") | str_detect(data$FL_187_DO, "D-R"), "Republican", NA)
data$party2_B <- ifelse(str_detect(data$FL_181_DO, "R-D") | str_detect(data$FL_187_DO, "R-D"), "Democrat", data$party2_B)
data$news2_B <- data$`G-2-2-1`
data$info_env2 <- ifelse(str_detect(data$FL_181_DO, "High") | str_detect(data$FL_187_DO, "High"), "High", NA)
data$info_env2 <- ifelse(str_detect(data$FL_181_DO, "Moderate") | str_detect(data$FL_187_DO, "Moderate"), "Moderate", data$info_env2)
data$info_env2 <- ifelse(str_detect(data$FL_181_DO, "Low") | str_detect(data$FL_187_DO, "Low"), "Low", data$info_env2)
data$choice2 <- ifelse(data$exp_arm=="Forced", paste0(data$Q112, data$Q116, data$Q120, data$Q124, data$Q128, data$Q132),
                       paste0(data$Q259, data$Q263, data$Q267, data$Q271, data$Q275, data$Q279))

data$party3_A <- ifelse(str_detect(data$FL_182_DO, "D-R") | str_detect(data$FL_188_DO, "D-R"), "Democrat", NA)
data$party3_A <- ifelse(str_detect(data$FL_182_DO, "R-D") | str_detect(data$FL_188_DO, "R-D"), "Republican", data$party3_A)
data$news3_A <- data$`G-3-1-1`
data$party3_B <- ifelse(str_detect(data$FL_182_DO, "D-R") | str_detect(data$FL_188_DO, "D-R"), "Republican", NA)
data$party3_B <- ifelse(str_detect(data$FL_182_DO, "R-D") | str_detect(data$FL_188_DO, "R-D"), "Democrat", data$party3_B)
data$news3_B <- data$`G-3-2-1`
data$info_env3 <- ifelse(str_detect(data$FL_182_DO, "High") | str_detect(data$FL_188_DO, "High"), "High", NA)
data$info_env3 <- ifelse(str_detect(data$FL_182_DO, "Moderate") | str_detect(data$FL_188_DO, "Moderate"), "Moderate", data$info_env3)
data$info_env3 <- ifelse(str_detect(data$FL_182_DO, "Low") | str_detect(data$FL_188_DO, "Low"), "Low", data$info_env3)
data$choice3 <- ifelse(data$exp_arm=="Forced", paste0(data$Q136, data$Q140, data$Q144, data$Q148, data$Q152, data$Q156),
                       paste0(data$Q283, data$Q287, data$Q291, data$Q295, data$Q299, data$Q303))

data$party4_A <- ifelse(str_detect(data$FL_183_DO, "D-R") | str_detect(data$FL_189_DO, "D-R"), "Democrat", NA)
data$party4_A <- ifelse(str_detect(data$FL_183_DO, "R-D") | str_detect(data$FL_189_DO, "R-D"), "Republican", data$party4_A)
data$news4_A <- data$`G-4-1-1`
data$party4_B <- ifelse(str_detect(data$FL_183_DO, "D-R") | str_detect(data$FL_189_DO, "D-R"), "Republican", NA)
data$party4_B <- ifelse(str_detect(data$FL_183_DO, "R-D") | str_detect(data$FL_189_DO, "R-D"), "Democrat", data$party4_B)
data$news4_B <- data$`G-4-2-1`
data$info_env4 <- ifelse(str_detect(data$FL_183_DO, "High") | str_detect(data$FL_189_DO, "High"), "High", NA)
data$info_env4 <- ifelse(str_detect(data$FL_183_DO, "Moderate") | str_detect(data$FL_189_DO, "Moderate"), "Moderate", data$info_env4)
data$info_env4 <- ifelse(str_detect(data$FL_183_DO, "Low") | str_detect(data$FL_189_DO, "Low"), "Low", data$info_env4)
data$choice4 <- ifelse(data$exp_arm=="Forced", paste0(data$Q162, data$Q166, data$Q170, data$Q178, data$Q174, data$Q182),
                       paste0(data$Q307, data$Q311, data$Q315, data$Q319, data$Q323, data$Q327))

data$party5_A <- ifelse(str_detect(data$FL_184_DO, "D-R") | str_detect(data$FL_190_DO, "D-R"), "Democrat", NA)
data$party5_A <- ifelse(str_detect(data$FL_184_DO, "R-D") | str_detect(data$FL_190_DO, "R-D"), "Republican", data$party5_A)
data$news5_A <- data$`G-5-1-1`
data$party5_B <- ifelse(str_detect(data$FL_184_DO, "D-R") | str_detect(data$FL_190_DO, "D-R"), "Republican", NA)
data$party5_B <- ifelse(str_detect(data$FL_184_DO, "R-D") | str_detect(data$FL_190_DO, "R-D"), "Democrat", data$party5_B)
data$news5_B <- data$`G-5-2-1`
data$info_env5 <- ifelse(str_detect(data$FL_184_DO, "High") | str_detect(data$FL_190_DO, "High"), "High", NA)
data$info_env5 <- ifelse(str_detect(data$FL_184_DO, "Moderate") | str_detect(data$FL_190_DO, "Moderate"), "Moderate", data$info_env5)
data$info_env5 <- ifelse(str_detect(data$FL_184_DO, "Low") | str_detect(data$FL_190_DO, "Low"), "Low", data$info_env5)
data$choice5 <- ifelse(data$exp_arm=="Forced", paste0(data$Q186, data$Q190, data$Q194, data$Q198, data$Q202, data$Q206),
                       paste0(data$Q331, data$Q335, data$Q339, data$Q343, data$Q347, data$Q351))

data$party6_A <- ifelse(str_detect(data$FL_185_DO, "D-R") | str_detect(data$FL_191_DO, "D-R"), "Democrat", NA)
data$party6_A <- ifelse(str_detect(data$FL_185_DO, "R-D") | str_detect(data$FL_191_DO, "R-D"), "Republican", data$party6_A)
data$news6_A <- data$`G-6-1-1`
data$party6_B <- ifelse(str_detect(data$FL_185_DO, "D-R") | str_detect(data$FL_191_DO, "D-R"), "Republican", NA)
data$party6_B <- ifelse(str_detect(data$FL_185_DO, "R-D") | str_detect(data$FL_191_DO, "R-D"), "Democrat", data$party6_B)
data$news6_B <- data$`G-6-2-1`
data$info_env6 <- ifelse(str_detect(data$FL_185_DO, "High") | str_detect(data$FL_191_DO, "High"), "High", NA)
data$info_env6 <- ifelse(str_detect(data$FL_185_DO, "Moderate") | str_detect(data$FL_191_DO, "Moderate"), "Moderate", data$info_env6)
data$info_env6 <- ifelse(str_detect(data$FL_185_DO, "Low") | str_detect(data$FL_191_DO, "Low"), "Low", data$info_env6)
data$choice6 <- ifelse(data$exp_arm=="Forced", paste0(data$Q210, data$Q214, data$Q218, data$Q226, data$Q222, data$Q230),
                       paste0(data$Q355, data$Q359, data$Q363, data$Q367, data$Q371, data$Q375))

# coding attribute-levels and choices for seventh task
data$party7_A <- ifelse(str_detect(data$FL_180_DO, "D-R") | str_detect(data$FL_186_DO, "D-R"), "Democrat", NA)
data$party7_A <- ifelse(str_detect(data$FL_180_DO, "R-D") | str_detect(data$FL_186_DO, "R-D"), "Republican", data$party7_A)
data$news7_A <- data$`G-1-1-1`
data$party7_B <- ifelse(str_detect(data$FL_180_DO, "D-R") | str_detect(data$FL_186_DO, "D-R"), "Republican", NA)
data$party7_B <- ifelse(str_detect(data$FL_180_DO, "R-D") | str_detect(data$FL_186_DO, "R-D"), "Democrat", data$party7_B)
data$news7_B <- data$`G-1-2-1`
data$info_env7 <- ifelse(str_detect(data$FL_180_DO, "High") | str_detect(data$FL_186_DO, "High"), "High", NA)
data$info_env7 <- ifelse(str_detect(data$FL_180_DO, "Moderate") | str_detect(data$FL_186_DO, "Moderate"), "Moderate", data$info_env7)
data$info_env7 <- ifelse(str_detect(data$FL_180_DO, "Low") | str_detect(data$FL_186_DO, "Low"), "Low", data$info_env7)
data$choice7 <- ifelse(data$exp_arm=="Abstention", paste0(data$Q86, data$Q92, data$Q96, data$Q104, data$Q100, data$Q108),
                       paste0(data$Q235, data$Q239, data$Q243, data$Q247, data$Q251, data$Q255))

# the following 12 code chunks reorganize this data from wide format to long format
# such that each row is a profile evaluated in a task by a respondent
long_data1_A <- subset(data, select = c(ResponseId, party1_A, news1_A, info_env1, choice1))
long_data1_A$task <- 1
long_data1_A$profile <- 1
colnames(long_data1_A) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data1_B <- subset(data, select = c(ResponseId, party1_B, news1_B, info_env1, choice1))
long_data1_B$task <- 1
long_data1_B$profile <- 2
colnames(long_data1_B) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data2_A <- subset(data, select = c(ResponseId, party2_A, news2_A, info_env2, choice2))
long_data2_A$task <- 2
long_data2_A$profile <- 1
colnames(long_data2_A) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data2_B <- subset(data, select = c(ResponseId, party2_B, news2_B, info_env2, choice2))
long_data2_B$task <- 2
long_data2_B$profile <- 2
colnames(long_data2_B) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data3_A <- subset(data, select = c(ResponseId, party3_A, news3_A, info_env3, choice3))
long_data3_A$task <- 3
long_data3_A$profile <- 1
colnames(long_data3_A) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data3_B <- subset(data, select = c(ResponseId, party3_B, news3_B, info_env3, choice3))
long_data3_B$task <- 3
long_data3_B$profile <- 2
colnames(long_data3_B) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data4_A <- subset(data, select = c(ResponseId, party4_A, news4_A, info_env4, choice4))
long_data4_A$task <- 4
long_data4_A$profile <- 1
colnames(long_data4_A) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data4_B <- subset(data, select = c(ResponseId, party4_B, news4_B, info_env4, choice4))
long_data4_B$task <- 4
long_data4_B$profile <- 2
colnames(long_data4_B) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data5_A <- subset(data, select = c(ResponseId, party5_A, news5_A, info_env5, choice5))
long_data5_A$task <- 5
long_data5_A$profile <- 1
colnames(long_data5_A) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data5_B <- subset(data, select = c(ResponseId, party5_B, news5_B, info_env5, choice5))
long_data5_B$task <- 5
long_data5_B$profile <- 2
colnames(long_data5_B) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data6_A <- subset(data, select = c(ResponseId, party6_A, news6_A, info_env6, choice6))
long_data6_A$task <- 6
long_data6_A$profile <- 1
colnames(long_data6_A) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data6_B <- subset(data, select = c(ResponseId, party6_B, news6_B, info_env6, choice6))
long_data6_B$task <- 6
long_data6_B$profile <- 2
colnames(long_data6_B) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data7_A <- subset(data, select = c(ResponseId, party7_A, news7_A, info_env7, choice7))
long_data7_A$task <- 7
long_data7_A$profile <- 1
colnames(long_data7_A) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

long_data7_B <- subset(data, select = c(ResponseId, party7_B, news7_B, info_env7, choice7))
long_data7_B$task <- 7
long_data7_B$profile <- 2
colnames(long_data7_B) <- c("ResponseId", "party", "news", "info_env", "choice", "task", "profile")

# this line binds the long format dataframes from each profile in each task
# into a single dataframe
long_data <- rbind(long_data1_A, long_data1_B, long_data2_A, long_data2_B, 
                   long_data3_A, long_data3_B, long_data4_A, long_data4_B,
                   long_data5_A, long_data5_B, long_data6_A, long_data6_B,
                   long_data7_A, long_data7_B)

# this line merges the respondent-level data back in to the long dataframe
full_data <- merge(demo_data, long_data, by = "ResponseId")

# the following five lines recode respondents' choices from character to numeric
full_data$choice_bin <- ifelse(full_data$choice=="Candidate A" & full_data$profile==1, 1, NA)
full_data$choice_bin <- ifelse(full_data$choice=="Candidate A" & full_data$profile==2, 0, full_data$choice_bin)
full_data$choice_bin <- ifelse(full_data$choice=="Candidate B" & full_data$profile==2, 1, full_data$choice_bin)
full_data$choice_bin <- ifelse(full_data$choice=="Candidate B" & full_data$profile==1, 0, full_data$choice_bin)
full_data$choice_bin <- ifelse(full_data$choice=="Neither candidate", 0, full_data$choice_bin)

# the following lines of code specify the levels and preferred ordering of the
# profile attribute-levels, information environment, and experimental arm
full_data$party <- factor(full_data$party, levels = c("Democrat", "Republican"))
full_data$news <- factor(full_data$news, levels = c("No recent news",
                                                    "Recently accused of sexual harassment",
                                                    "Recently accused of leaking confidential information",
                                                    "Recently accused of cheating on spouse",
                                                    "Recently celebrated wedding anniversary",
                                                    "Recently honored for public service"))
full_data$info_env <- factor(full_data$info_env, levels = c("Low",
                                                            "Moderate",
                                                            "High"))
full_data$exp_arm <- factor(full_data$exp_arm, levels = c("Forced", "Abstention"))

# The following lines estimate our primary model of interest, which interacts
# respondents' experimental arm, their information environment, and the recent
# news attribute-level for the profile of interest and then uses the
# cluster.vcov function from the multiwayvcov package to obtain the variance-
# covariance matrix that corresponds with clustering at the respondent-level.
# We limit ourselves to tasks 1-6, as respondents switched experimental arms for
# the seventh task

# NOTE: this object is titled "table2_arm" because it is modeled after Funck and
# McCabe's Table 2
 
table2_arm <- lm(choice_bin ~ info_env*news*exp_arm, data = full_data[which(full_data$task!=7)])
table2_arm_cluster_vcov <- cluster.vcov(table2_arm, cluster=full_data$ResponseId[which(full_data$task!=7)])
coeftest(table2_arm, table2_arm_cluster_vcov)
table2_arm_ses <- coeftest(table2_arm, table2_arm_cluster_vcov)[,2]
table2_arm_pvals <- coeftest(table2_arm, table2_arm_cluster_vcov)[,4]

# This line of code replaces standard variance-covariance matrix with the
# cluster robust variance-covariance
table2_arm <-withVCov(table2_arm, vcov = table2_arm_cluster_vcov)

# In the below lines of code, we calculate the effect of the non-baseline levels
# of the recent news attribute in the moderate and high information environments
# relative to that for the low information environment.  The first 10 estimates
# correspond with the forced choice arm, while the next 10 estimates correspond
# with the abstention arm.  In accordance with our pre-analysis plan, we
# use a Bonferroni-corrected critical value of 2.81.
diffs_from_low_info <- confint(glht(table2_arm,
                                    linfct = c("`info_envModerate:newsRecently accused of sexual harassment` = 0",
                                               "`info_envModerate:newsRecently accused of leaking confidential information` = 0",
                                               "`info_envModerate:newsRecently accused of cheating on spouse` = 0",
                                               "`info_envModerate:newsRecently celebrated wedding anniversary` = 0",
                                               "`info_envModerate:newsRecently honored for public service` = 0",
                                               "`info_envHigh:newsRecently accused of sexual harassment` = 0",
                                               "`info_envHigh:newsRecently accused of leaking confidential information` = 0",
                                               "`info_envHigh:newsRecently accused of cheating on spouse` = 0",
                                               "`info_envHigh:newsRecently celebrated wedding anniversary` = 0",
                                               "`info_envHigh:newsRecently honored for public service` = 0",
                                               "`info_envModerate:newsRecently accused of sexual harassment` + `info_envModerate:newsRecently accused of sexual harassment:exp_armAbstention` = 0",
                                               "`info_envModerate:newsRecently accused of leaking confidential information` + `info_envModerate:newsRecently accused of leaking confidential information:exp_armAbstention` = 0",
                                               "`info_envModerate:newsRecently accused of cheating on spouse` + `info_envModerate:newsRecently accused of cheating on spouse:exp_armAbstention` = 0",
                                               "`info_envModerate:newsRecently celebrated wedding anniversary` + `info_envModerate:newsRecently celebrated wedding anniversary:exp_armAbstention` = 0",
                                               "`info_envModerate:newsRecently honored for public service` + `info_envModerate:newsRecently honored for public service:exp_armAbstention` = 0",
                                               "`info_envHigh:newsRecently accused of sexual harassment` + `info_envHigh:newsRecently accused of sexual harassment:exp_armAbstention` = 0",
                                               "`info_envHigh:newsRecently accused of leaking confidential information` + `info_envHigh:newsRecently accused of leaking confidential information:exp_armAbstention` = 0",
                                               "`info_envHigh:newsRecently accused of cheating on spouse` + `info_envHigh:newsRecently accused of cheating on spouse:exp_armAbstention` = 0",
                                               "`info_envHigh:newsRecently celebrated wedding anniversary` + `info_envHigh:newsRecently celebrated wedding anniversary:exp_armAbstention` = 0",
                                               "`info_envHigh:newsRecently honored for public service` + `info_envHigh:newsRecently honored for public service:exp_armAbstention` = 0")),
                               calpha = 2.81)$confint # 20 hypothesis tests, 0.05/20=0.0025, corresponding critical value 2.81

# In the below lines of code, we calculate the differences between the estimates
# obtained for each recent news attribute-level in the low vs. moderate and
# high information environments across the abstention and forced choice arms
diffs_in_diffs <- confint(glht(table2_arm,
                               linfct = c("`info_envModerate:newsRecently accused of sexual harassment:exp_armAbstention` = 0",
                                          "`info_envModerate:newsRecently accused of leaking confidential information:exp_armAbstention` = 0",
                                          "`info_envModerate:newsRecently accused of cheating on spouse:exp_armAbstention` = 0",
                                          "`info_envModerate:newsRecently celebrated wedding anniversary:exp_armAbstention` = 0",
                                          "`info_envModerate:newsRecently honored for public service:exp_armAbstention` = 0",
                                          "`info_envHigh:newsRecently accused of sexual harassment:exp_armAbstention` = 0",
                                          "`info_envHigh:newsRecently accused of leaking confidential information:exp_armAbstention` = 0",
                                          "`info_envHigh:newsRecently accused of cheating on spouse:exp_armAbstention` = 0",
                                          "`info_envHigh:newsRecently celebrated wedding anniversary:exp_armAbstention` = 0",
                                          "`info_envHigh:newsRecently honored for public service:exp_armAbstention` = 0")),
                          calpha = 2.81)$confint # 20 hypothesis tests, 0.05/20=0.0025, corresponding critical value 2.81

# we multiply the differences in the differences by -1 to reflect that we want to
# express the following quantity: forced choice AMCE - abstention AMCE
diffs_in_diffs <- -1*diffs_in_diffs

##########
# FIGURE 1
##########

pdf(file = "../figures/figure1.pdf", family = "Times", height = 11, width=9)
layout(matrix(c(1,2,3,3), nrow=2, byrow=TRUE), heights = c(0.90,0.10),
       widths = c(0.6,0.4))
par(mar=c(6.1,12,1.75,1))
plot(x=c(diffs_from_low_info[11,1],
         diffs_from_low_info[1,1],
         diffs_from_low_info[12,1],
         diffs_from_low_info[2,1],
         diffs_from_low_info[13,1],
         diffs_from_low_info[3,1],
         diffs_from_low_info[14,1],
         diffs_from_low_info[4,1],
         diffs_from_low_info[15,1],
         diffs_from_low_info[5,1],
         diffs_from_low_info[16,1],
         diffs_from_low_info[6,1],
         diffs_from_low_info[17,1],
         diffs_from_low_info[7,1],
         diffs_from_low_info[18,1],
         diffs_from_low_info[8,1],
         diffs_from_low_info[19,1],
         diffs_from_low_info[9,1],
         diffs_from_low_info[20,1],
         diffs_from_low_info[10,1]), 
     y=c(1.3, 1.7,
         2.3, 2.7,
         3.3, 3.7,
         4.3, 4.7,
         5.3, 5.7,
         7.3, 7.7,
         8.3, 8.7,
         9.3, 9.7,
         10.3, 10.7,
         11.3, 11.7),
     xlim=c(-0.25,0.25),
     xaxt="n",
     ylim=c(1, 12),
     pch=c(rep(c(17,19))),
     col=c(rep(c("gray50","black"))),
     
     tck=-.02,
     cex.axis=0.9,
     cex=1.25,
     ylab="",
     yaxt="n",
     xlab="",
     axes = FALSE,
     panel.first = c(abline(v=0,lwd=2, col="gray70",lty=2)))
segments(col=c(rep(c("gray50","black"),18)),
         x0=c(diffs_from_low_info[11,2],
              diffs_from_low_info[1,2],
              diffs_from_low_info[12,2],
              diffs_from_low_info[2,2],
              diffs_from_low_info[13,2],
              diffs_from_low_info[3,2],
              diffs_from_low_info[14,2],
              diffs_from_low_info[4,2],
              diffs_from_low_info[15,2],
              diffs_from_low_info[5,2],
              diffs_from_low_info[16,2],
              diffs_from_low_info[6,2],
              diffs_from_low_info[17,2],
              diffs_from_low_info[7,2],
              diffs_from_low_info[18,2],
              diffs_from_low_info[8,2],
              diffs_from_low_info[19,2],
              diffs_from_low_info[9,2],
              diffs_from_low_info[20,2],
              diffs_from_low_info[10,2]),
         x1=c(diffs_from_low_info[11,3],
              diffs_from_low_info[1,3],
              diffs_from_low_info[12,3],
              diffs_from_low_info[2,3],
              diffs_from_low_info[13,3],
              diffs_from_low_info[3,3],
              diffs_from_low_info[14,3],
              diffs_from_low_info[4,3],
              diffs_from_low_info[15,3],
              diffs_from_low_info[5,3],
              diffs_from_low_info[16,3],
              diffs_from_low_info[6,3],
              diffs_from_low_info[17,3],
              diffs_from_low_info[7,3],
              diffs_from_low_info[18,3],
              diffs_from_low_info[8,3],
              diffs_from_low_info[19,3],
              diffs_from_low_info[9,3],
              diffs_from_low_info[20,3],
              diffs_from_low_info[10,3]),
         y0=c(1.3, 1.7,
              2.3, 2.7,
              3.3, 3.7,
              4.3, 4.7,
              5.3, 5.7,
              7.3, 7.7,
              8.3, 8.7,
              9.3, 9.7,
              10.3, 10.7,
              11.3, 11.7
         ))
axis(1,at=c(seq(-0.25,0.25,0.05)), 
     labels = c("-25%", NA, NA, NA, NA, "0%", NA, NA, NA, NA, "25%"), cex.axis = 1.5)
par(mgp=c(0,8,0))
axis(2,at=c(seq(1.5,5.5,1),
            seq(7.5,11.5,1)),
     las=2,
     labels=rep(c("Sexual Harassment",
                  "Leaking Confidential\nInformation",
                  "Cheating on Spouse",
                  "Celebrated Wedding\nAnniversary",
                  "Honored for Public\nService"),2),
     tck=0,
     lwd = 0,
     line = 1,
     cex.axis=1, hadj=0)
brackets(x1=-.26, x2=-.26, y1=1.1, y2=5.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Moderate Information")), side = 2, at=3.5, line = 10.2, cex = 1.75)
brackets(x1=-.26, x2=-.26, y1=7.1, y2=11.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("High Information")), side = 2, at=9.5, line = 10.2, cex = 1.75)
par(mgp=c(0,13,0))
mtext("Difference from\nLow Information AMCE", side = 1, at=0, line = 5, cex = 1.75)


par(mar=c(6.1,1,1.75,1),
    mgp=c(3,1,0))
plot(x=c(diffs_in_diffs[1,1],
         diffs_in_diffs[2,1],
         diffs_in_diffs[3,1],
         diffs_in_diffs[4,1],
         diffs_in_diffs[5,1],
         diffs_in_diffs[6,1],
         diffs_in_diffs[7,1],
         diffs_in_diffs[8,1],
         diffs_in_diffs[9,1],
         diffs_in_diffs[10,1]), 
     y=c(1.5,
         2.5,
         3.5,
         4.5,
         5.5,
         7.5,
         8.5,
         9.5,
         10.5,
         11.5),
     xlim=c(-0.25,0.25),
     xaxt="n",
     ylim=c(1, 12),
     pch=15,
     col="black",
     tck=-.02,
     cex.axis=0.9,
     cex=1.25,
     ylab="",
     yaxt="n",
     xlab="",
     axes = FALSE,
     panel.first = c(abline(v=0,lwd=2, col="gray70",lty=2)))
segments(col="black",
         x0=c(diffs_in_diffs[1,2],
              diffs_in_diffs[2,2],
              diffs_in_diffs[3,2],
              diffs_in_diffs[4,2],
              diffs_in_diffs[5,2],
              diffs_in_diffs[6,2],
              diffs_in_diffs[7,2],
              diffs_in_diffs[8,2],
              diffs_in_diffs[9,2],
              diffs_in_diffs[10,2]),
         x1=c(diffs_in_diffs[1,3],
              diffs_in_diffs[2,3],
              diffs_in_diffs[3,3],
              diffs_in_diffs[4,3],
              diffs_in_diffs[5,3],
              diffs_in_diffs[6,3],
              diffs_in_diffs[7,3],
              diffs_in_diffs[8,3],
              diffs_in_diffs[9,3],
              diffs_in_diffs[10,3]),
         y0=c(1.5,
              2.5,
              3.5,
              4.5,
              5.5,
              7.5,
              8.5,
              9.5,
              10.5,
              11.5
         ))
axis(1,at=c(seq(-0.25,0.25,0.05)), 
     labels = c("-25%", NA, NA, NA, NA, "0%", NA, NA, NA, NA, "25%"), cex.axis = 1.5)
par(mgp=c(0,13,0))
mtext("Forced Choice AMCE -\nAbstention Option AMCE", side = 1, at=0, line = 5, cex = 1.75)

par(mar=c(0,0,0,0))
plot(0,0, type="n", axes=FALSE, xlab="", ylab="")
legend(x=-0.06, y=0.8, legend=c("Forced Choice",
                                "Abstention Option"), lty=1, lwd=1, 
       pch = c(19,17),
       col = c("black", "grey50"), cex = 1.50)
dev.off()

############
# TABLE SM.5
############

# here, we mirror Funck and McCabe's presentation of results in their main
# Table 1 by separately estimating a model where the reference category is
# low information (vs. moderate or high) and is moderate information (vs. high).
# We then additional present these results separately for our abstention and
# forcec choice arms.

# in each chunk of code, we estimate the model with the appropriate baseline
# category, then obtain the cluster robust variance-covariance matrix (clustered
# on respondent) and save the corresponding standard errors and p-values

table2_forced <- lm(choice_bin ~ info_env*news, data = full_data[which(full_data$exp_arm=="Forced" & full_data$task!=7),])
table2_forced_cluster_vcov <- cluster.vcov(table2_forced, cluster=full_data$ResponseId[which(full_data$exp_arm=="Forced" & full_data$task!=7)])
coeftest(table2_forced, table2_forced_cluster_vcov)
table2_forced_ses <- coeftest(table2_forced, table2_forced_cluster_vcov)[,2]
table2_forced_pvals <- coeftest(table2_forced, table2_forced_cluster_vcov)[,4]

table2_forced_medium <- lm(choice_bin ~ info_env*news, data = full_data[which(full_data$exp_arm=="Forced" & full_data$task!=7 & full_data$info_env!="Low"),])
table2_forced_medium_cluster_vcov <- cluster.vcov(table2_forced_medium, cluster=full_data$ResponseId[which(full_data$exp_arm=="Forced" & full_data$task!=7 & full_data$info_env!="Low")])
coeftest(table2_forced_medium, table2_forced_medium_cluster_vcov)
table2_forced_medium_ses <- coeftest(table2_forced_medium, table2_forced_medium_cluster_vcov)[,2]
table2_forced_medium_pvals <- coeftest(table2_forced_medium, table2_forced_medium_cluster_vcov)[,4]

table2_abstention <- lm(choice_bin ~ info_env*news, data = full_data[which(full_data$exp_arm=="Abstention" & full_data$task!=7),])
table2_abstention_cluster_vcov <- cluster.vcov(table2_abstention, cluster=full_data$ResponseId[which(full_data$exp_arm=="Abstention" & full_data$task!=7)])
coeftest(table2_abstention, table2_abstention_cluster_vcov)
table2_abstention_ses <- coeftest(table2_abstention, table2_abstention_cluster_vcov)[,2]
table2_abstention_pvals <- coeftest(table2_abstention, table2_abstention_cluster_vcov)[,4]

table2_abstention_medium <- lm(choice_bin ~ info_env*news, data = full_data[which(full_data$exp_arm=="Abstention" & full_data$task!=7 & full_data$info_env!="Low"),])
table2_abstention_medium_cluster_vcov <- cluster.vcov(table2_abstention_medium, cluster=full_data$ResponseId[which(full_data$exp_arm=="Abstention" & full_data$task!=7 & full_data$info_env!="Low")])
coeftest(table2_abstention_medium, table2_abstention_medium_cluster_vcov)
table2_abstention_medium_ses <- coeftest(table2_abstention_medium, table2_abstention_medium_cluster_vcov)[,2]
table2_abstention_medium_pvals <- coeftest(table2_abstention_medium, table2_abstention_medium_cluster_vcov)[,4]

print(
  texreg(l=list(table2_forced, table2_forced_medium,
              table2_abstention, table2_abstention_medium),
       custom.coef.map = list("(Intercept)" = "Intercept",
                              "info_envModerate" = "Medium information",
                              "info_envHigh" = "High information",
                              "newsRecently accused of cheating on spouse" = "Cheating on spouse",
                              "newsRecently accused of leaking confidential information" = "Leaking confidential information",
                              "newsRecently accused of sexual harassment" = "Sexual harassment",
                              "newsRecently celebrated wedding anniversary" = "Wedding anniversary",
                              "newsRecently honored for public service" = "Honored public service",
                              "info_envModerate:newsRecently accused of cheating on spouse" = "Medium information x cheating on spouse",
                              "info_envModerate:newsRecently accused of leaking confidential information" = "Medium information x leaking confidential information",
                              "info_envModerate:newsRecently accused of sexual harassment" = "Medium information x sexual harassment",
                              "info_envModerate:newsRecently celebrated wedding anniversary" = "Medium information x wedding anniversary",
                              "info_envModerate:newsRecently honored for public service" = "Medium information x honored public service",
                              "info_envHigh:newsRecently accused of cheating on spouse" = "High information x cheating on spouse",
                              "info_envHigh:newsRecently accused of leaking confidential information" = "High information x leaking confidential information",
                              "info_envHigh:newsRecently accused of sexual harassment" = "High information x sexual harassment",
                              "info_envHigh:newsRecently celebrated wedding anniversary" = "High information x wedding anniversary",
                              "info_envHigh:newsRecently honored for public service" = "High information x honored public service"),
       custom.model.names = c("Ref. Cat. = Low", "Ref. Cat. = Medium", 
                              "Ref. Cat. = Low", "Ref. Cat. = Medium"),
       override.se = list(table2_forced_ses, table2_forced_medium_ses,
                          table2_abstention_ses, table2_abstention_medium_ses),
       override.pvalues = list(table2_forced_pvals, table2_forced_medium_pvals,
                               table2_abstention_pvals, table2_abstention_medium_pvals),
       caption = "Funck and McCabe (2021) Replication Results (Table 2)",
       caption.above = TRUE,
       stars = 0.05,
       include.rsquared = FALSE,
       include.adjrs = FALSE)
  , file = "../tables/table_SM5.tex"
)

#############
# FIGURE SM.2
#############

# calculating difference in AMCE between low information and moderate and high
# information environments
info_env_amces <- confint(glht(table2_arm, 
                               linfct = c("`newsRecently accused of sexual harassment` = 0",
                                          "`newsRecently accused of leaking confidential information` = 0",
                                          "`newsRecently accused of cheating on spouse` = 0",
                                          "`newsRecently celebrated wedding anniversary` = 0",
                                          "`newsRecently honored for public service` = 0",
                                          "`newsRecently accused of sexual harassment` + `info_envModerate:newsRecently accused of sexual harassment` = 0",
                                          "`newsRecently accused of leaking confidential information` + `info_envModerate:newsRecently accused of leaking confidential information`= 0",
                                          "`newsRecently accused of cheating on spouse` + `info_envModerate:newsRecently accused of cheating on spouse` = 0",
                                          "`newsRecently celebrated wedding anniversary` + `info_envModerate:newsRecently celebrated wedding anniversary` = 0",
                                          "`newsRecently honored for public service` + `info_envModerate:newsRecently honored for public service` = 0",
                                          "`newsRecently accused of sexual harassment` + `info_envHigh:newsRecently accused of sexual harassment` = 0",
                                          "`newsRecently accused of leaking confidential information` + `info_envHigh:newsRecently accused of leaking confidential information`= 0",
                                          "`newsRecently accused of cheating on spouse` + `info_envHigh:newsRecently accused of cheating on spouse` = 0",
                                          "`newsRecently celebrated wedding anniversary` + `info_envHigh:newsRecently celebrated wedding anniversary` = 0",
                                          "`newsRecently honored for public service` + `info_envHigh:newsRecently honored for public service` = 0",
                                          "`newsRecently accused of sexual harassment` + `newsRecently accused of sexual harassment:exp_armAbstention` = 0",
                                          "`newsRecently accused of leaking confidential information` + `newsRecently accused of leaking confidential information:exp_armAbstention` = 0",
                                          "`newsRecently accused of cheating on spouse` + `newsRecently accused of cheating on spouse:exp_armAbstention` = 0",
                                          "`newsRecently celebrated wedding anniversary` + `newsRecently celebrated wedding anniversary:exp_armAbstention` = 0",
                                          "`newsRecently honored for public service` + `newsRecently honored for public service:exp_armAbstention` = 0",
                                          "`newsRecently accused of sexual harassment` + `info_envModerate:newsRecently accused of sexual harassment` + `newsRecently accused of sexual harassment:exp_armAbstention` + `info_envModerate:newsRecently accused of sexual harassment:exp_armAbstention` = 0",
                                          "`newsRecently accused of leaking confidential information` + `info_envModerate:newsRecently accused of leaking confidential information` + `newsRecently accused of leaking confidential information:exp_armAbstention` + `info_envModerate:newsRecently accused of leaking confidential information:exp_armAbstention` = 0",
                                          "`newsRecently accused of cheating on spouse` + `info_envModerate:newsRecently accused of cheating on spouse` + `newsRecently accused of cheating on spouse:exp_armAbstention` + `info_envModerate:newsRecently accused of cheating on spouse:exp_armAbstention` = 0",
                                          "`newsRecently celebrated wedding anniversary` + `info_envModerate:newsRecently celebrated wedding anniversary` + `newsRecently celebrated wedding anniversary:exp_armAbstention` + `info_envModerate:newsRecently celebrated wedding anniversary:exp_armAbstention` = 0",
                                          "`newsRecently honored for public service` + `info_envModerate:newsRecently honored for public service` + `newsRecently honored for public service:exp_armAbstention` + `info_envModerate:newsRecently honored for public service:exp_armAbstention` = 0",
                                          "`newsRecently accused of sexual harassment` + `info_envHigh:newsRecently accused of sexual harassment` +`newsRecently accused of sexual harassment:exp_armAbstention` + `info_envHigh:newsRecently accused of sexual harassment:exp_armAbstention` = 0",
                                          "`newsRecently accused of leaking confidential information` + `info_envHigh:newsRecently accused of leaking confidential information` + `newsRecently accused of leaking confidential information:exp_armAbstention` + `info_envHigh:newsRecently accused of leaking confidential information:exp_armAbstention` = 0",
                                          "`newsRecently accused of cheating on spouse` + `info_envHigh:newsRecently accused of cheating on spouse` + `newsRecently accused of cheating on spouse:exp_armAbstention` + `info_envHigh:newsRecently accused of cheating on spouse:exp_armAbstention` = 0",
                                          "`newsRecently celebrated wedding anniversary` + `info_envHigh:newsRecently celebrated wedding anniversary` + `newsRecently celebrated wedding anniversary:exp_armAbstention` + `info_envHigh:newsRecently celebrated wedding anniversary:exp_armAbstention` = 0",
                                          "`newsRecently honored for public service` + `info_envHigh:newsRecently honored for public service` + `newsRecently honored for public service:exp_armAbstention` + `info_envHigh:newsRecently honored for public service:exp_armAbstention` = 0")),
                          calpha = univariate_calpha())$confint

# calculating difference in differences between low and moderate and high
# information environments across abstention and forced choice arms
diffs_by_info_env <- confint(glht(table2_arm, 
                                  linfct = c("`newsRecently accused of sexual harassment:exp_armAbstention` = 0",
                                             "`newsRecently accused of leaking confidential information:exp_armAbstention` = 0",
                                             "`newsRecently accused of cheating on spouse:exp_armAbstention` = 0",
                                             "`newsRecently celebrated wedding anniversary:exp_armAbstention` = 0",
                                             "`newsRecently honored for public service:exp_armAbstention` = 0",
                                             "`newsRecently accused of sexual harassment:exp_armAbstention` + `info_envModerate:newsRecently accused of sexual harassment:exp_armAbstention` = 0",
                                             "`newsRecently accused of leaking confidential information:exp_armAbstention` + `info_envModerate:newsRecently accused of leaking confidential information:exp_armAbstention` = 0",
                                             "`newsRecently accused of cheating on spouse:exp_armAbstention` + `info_envModerate:newsRecently accused of cheating on spouse:exp_armAbstention` = 0",
                                             "`newsRecently celebrated wedding anniversary:exp_armAbstention` + `info_envModerate:newsRecently celebrated wedding anniversary:exp_armAbstention` = 0",
                                             "`newsRecently honored for public service:exp_armAbstention` + `info_envModerate:newsRecently honored for public service:exp_armAbstention` = 0",
                                             "`newsRecently accused of sexual harassment:exp_armAbstention` + `info_envHigh:newsRecently accused of sexual harassment:exp_armAbstention` = 0",
                                             "`newsRecently accused of leaking confidential information:exp_armAbstention` + `info_envHigh:newsRecently accused of leaking confidential information:exp_armAbstention` = 0",
                                             "`newsRecently accused of cheating on spouse:exp_armAbstention` + `info_envHigh:newsRecently accused of cheating on spouse:exp_armAbstention` = 0",
                                             "`newsRecently celebrated wedding anniversary:exp_armAbstention` + `info_envHigh:newsRecently celebrated wedding anniversary:exp_armAbstention` = 0",
                                             "`newsRecently honored for public service:exp_armAbstention` + `info_envHigh:newsRecently honored for public service:exp_armAbstention` = 0")),
                             calpha = 2.71)$confint

# changing the sign of a few of the calculated distances to better represent the
# quantity we want to display--whether the magnitude of the forced choice estimate
# is larger (positive) or smaller (negative)
diffs_by_info_env[4,] <- -diffs_by_info_env[4,]
diffs_by_info_env[5,] <- -diffs_by_info_env[5,]
diffs_by_info_env[9,] <- -diffs_by_info_env[9,]
diffs_by_info_env[10,] <- -diffs_by_info_env[10,]

# reproducing Figure SM.2
pdf(file = "../figures/figure_SM2.pdf", family = "Times", height = 18, width=13)
layout(matrix(c(1,2), nrow=2, byrow=TRUE), heights = c(0.85,0.15))
par(mar=c(6.1,12,5.75,1))
plot(x=c(# low info
    info_env_amces[16,1],
    info_env_amces[1,1],
    info_env_amces[17,1],
    info_env_amces[2,1],
    info_env_amces[18,1],
    info_env_amces[3,1],
    0,
    0,
    info_env_amces[19,1],
    info_env_amces[4,1],
    info_env_amces[20,1],
    info_env_amces[5,1],
    # mod info
    info_env_amces[21,1],
    info_env_amces[6,1],
    info_env_amces[22,1],
    info_env_amces[7,1],
    info_env_amces[23,1],
    info_env_amces[8,1],
    0,
    0,
    info_env_amces[24,1],
    info_env_amces[9,1],
    info_env_amces[25,1],
    info_env_amces[10,1],
    # high info
    info_env_amces[26,1],
    info_env_amces[11,1],
    info_env_amces[27,1],
    info_env_amces[12,1],
    info_env_amces[28,1],
    info_env_amces[13,1],
    0,
    0,
    info_env_amces[29,1],
    info_env_amces[14,1],
    info_env_amces[30,1],
    info_env_amces[15,1],
    diffs_by_info_env[1:15,1]+1.15), 
    y=c(1.3, 1.7,
        2.3, 2.7,
        3.3, 3.7,
        4.3, 4.7,
        5.3, 5.7,
        6.3, 6.7,
        8.3, 8.7,
        9.3, 9.7,
        10.3, 10.7,
        11.3, 11.7,
        12.3, 12.7,
        13.3, 13.7,
        15.3, 15.7,
        16.3, 16.7,
        17.3, 17.7,
        18.3, 18.7,
        19.3, 19.7,
        20.3, 20.7,
        seq(1.5,3.5,1),
        seq(5.5,6.5,1),
        seq(8.5,10.5,1),
        seq(12.5,13.5,1),
        seq(15.5,17.5,1),
        seq(19.5,20.5,1)
    ),
    xlim=c(-0.5,1.40),
    xaxt="n",
    ylim=c(1, 21),
    pch=c(rep(c(17,19),18), rep(15, 16)),
    col=c(rep(c("gray50","black"),18), rep("black", 18)),
    
    tck=-.02,
    cex.axis=0.9,
    cex=1.25,
    ylab="",
    yaxt="n",
    xlab="",
    axes = FALSE,
    panel.first = c(abline(v=0,lwd=2, col="gray70",lty=2),
                    abline(v=1.15,lwd=2, col="gray70",lty=2)))
segments(col=c(rep(c("gray50","black"),18), rep("black", 18)),
         x0=c(# low info
             info_env_amces[16,2],
             info_env_amces[1,2],
             info_env_amces[17,2],
             info_env_amces[2,2],
             info_env_amces[18,2],
             info_env_amces[3,2],
             0,
             0,
             info_env_amces[19,2],
             info_env_amces[4,2],
             info_env_amces[20,2],
             info_env_amces[5,2],
             # mod info
             info_env_amces[21,2],
             info_env_amces[6,2],
             info_env_amces[22,2],
             info_env_amces[7,2],
             info_env_amces[23,2],
             info_env_amces[8,2],
             0,
             0,
             info_env_amces[24,2],
             info_env_amces[9,2],
             info_env_amces[25,2],
             info_env_amces[10,2],
             # high info
             info_env_amces[26,2],
             info_env_amces[11,2],
             info_env_amces[27,2],
             info_env_amces[12,2],
             info_env_amces[28,2],
             info_env_amces[13,2],
             0,
             0,
             info_env_amces[29,2],
             info_env_amces[14,2],
             info_env_amces[30,2],
             info_env_amces[15,2],
             diffs_by_info_env[1:15,2]+1.15
         ),
         x1=c(# low info
             info_env_amces[16,3],
             info_env_amces[1,3],
             info_env_amces[17,3],
             info_env_amces[2,3],
             info_env_amces[18,3],
             info_env_amces[3,3],
             0,
             0,
             info_env_amces[19,3],
             info_env_amces[4,3],
             info_env_amces[20,3],
             info_env_amces[5,3],
             # mod info
             info_env_amces[21,3],
             info_env_amces[6,3],
             info_env_amces[22,3],
             info_env_amces[7,3],
             info_env_amces[23,3],
             info_env_amces[8,3],
             0,
             0,
             info_env_amces[24,3],
             info_env_amces[9,3],
             info_env_amces[25,3],
             info_env_amces[10,3],
             # high info
             info_env_amces[26,3],
             info_env_amces[11,3],
             info_env_amces[27,3],
             info_env_amces[12,3],
             info_env_amces[28,3],
             info_env_amces[13,3],
             0,
             0,
             info_env_amces[29,3],
             info_env_amces[14,3],
             info_env_amces[30,3],
             info_env_amces[15,3],
             diffs_by_info_env[1:15,3]+1.15),
         y0=c(1.3, 1.7,
              2.3, 2.7,
              3.3, 3.7,
              4.3, 4.7,
              5.3, 5.7,
              6.3, 6.7,
              8.3, 8.7,
              9.3, 9.7,
              10.3, 10.7,
              11.3, 11.7,
              12.3, 12.7,
              13.3, 13.7,
              15.3, 15.7,
              16.3, 16.7,
              17.3, 17.7,
              18.3, 18.7,
              19.3, 19.7,
              20.3, 20.7,
              seq(1.5,3.5,1),
              seq(5.5,6.5,1),
              seq(8.5,10.5,1),
              seq(12.5,13.5,1),
              seq(15.5,17.5,1),
              seq(19.5,20.5,1)
         ))
axis(1,at=c(-0.50, -.25, 0, .25, .5), 
     labels = c("-50%", "-25%", "0%", "25%", "50%"), cex.axis = 1.5)
axis(1,at=c(.90, 1.15, 1.40), 
     labels = c("-25%", "0%", "25%"), cex.axis = 1.5)
par(mgp=c(0,8,0))
axis(2,at=c(seq(1.5,6.5,1),
            seq(8.5,13.5,1),
            seq(15.5,20.5,1)),
     las=2,
     labels=rep(c("Sexual Harassment",
                  "Leaking Confidential\nInformation",
                  "Cheating on Spouse",
                  "No Recent News",
                  "Celebrated Wedding\nAnniversary",
                  "Honored for Public\nService"),3),
     tck=0,
     lwd = 0,
     line = 0,
     cex.axis=1, hadj=0)
brackets(x1=-.53, x2=-.53, y1=1.1, y2=6.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Low Information")), side = 2, at=4, line = 9.2, cex = 1.75)
brackets(x1=-.53, x2=-.53, y1=8.1, y2=13.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Moderate Information")), side = 2, at=11, line = 9.2, cex = 1.75)
brackets(x1=-.53, x2=-.53, y1=15.1, y2=20.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("High Information")), side = 2, at=18, line = 9.2, cex = 1.75)
par(mgp=c(0,11,0))
mtext("Diff. from Baseline", side = 1, at=0, line = 4, cex = 1.75)
mtext("Magnitude of\nUnder/Overestimation", side = 1, at=1.15, line = 4, cex = 1.75)
mtext("AMCEs by\n Outcome Options", side = 3, at=0, line=0.5, cex = 2)
mtext("Under/Overestimation\nwith Forced Choice", side = 3, at=1.15, line=0.5, cex = 2)
par(mar=c(0,0,0,0))
plot(0,0, type="n", axes=FALSE, xlab="", ylab="")
legend(x=0, y=0.8, legend=c("Forced Choice",
                            "Abstention Option",
                            "Difference in AMCEs"), lty=1, lwd=1, 
       pch = c(19,17,15),
       col = c("black", "grey50", "black"), cex = 1.25)
dev.off()

############
# TABLE SM.6
############

# recoding respondent choices that were "Neither candidate" to be NA instead of
# 0
full_data$choice_bin_NA <- ifelse(full_data$choice=="Neither candidate", NA, full_data$choice_bin)

# as with our code for our Table SM.5, we mirror Funck and McCabe's presentation 
# of results in their main Table 1 by separately estimating a model where the 
# reference category is low information (vs. moderate or high) and is moderate 
# information (vs. high). We then additional present these results separately 
# for our abstention and forced choice arms.

# in each chunk of code, we estimate the model with the appropriate baseline
# category, then obtain the cluster robust variance-covariance matrix (clustered
# on respondent) and save the corresponding standard errors and p-values

table2_forced_NA <- lm(choice_bin_NA ~ info_env*news, data = full_data[which(full_data$exp_arm=="Forced" & full_data$task!=7),])
table2_forced_NA_cluster_vcov <- cluster.vcov(table2_forced_NA, cluster=full_data$ResponseId[which(full_data$exp_arm=="Forced" & full_data$task!=7)])
coeftest(table2_forced_NA, table2_forced_NA_cluster_vcov)
table2_forced_NA_ses <- coeftest(table2_forced_NA, table2_forced_NA_cluster_vcov)[,2]
table2_forced_NA_pvals <- coeftest(table2_forced_NA, table2_forced_NA_cluster_vcov)[,4]

table2_forced_NA_medium <- lm(choice_bin_NA ~ info_env*news, data = full_data[which(full_data$exp_arm=="Forced" & full_data$task!=7 & full_data$info_env!="Low"),])
table2_forced_NA_medium_cluster_vcov <- cluster.vcov(table2_forced_NA_medium, cluster=full_data$ResponseId[which(full_data$exp_arm=="Forced" & full_data$task!=7 & full_data$info_env!="Low")])
coeftest(table2_forced_NA_medium, table2_forced_NA_medium_cluster_vcov)
table2_forced_NA_medium_ses <- coeftest(table2_forced_NA_medium, table2_forced_NA_medium_cluster_vcov)[,2]
table2_forced_NA_medium_pvals <- coeftest(table2_forced_NA_medium, table2_forced_NA_medium_cluster_vcov)[,4]

table2_abstention_NA <- lm(choice_bin_NA ~ info_env*news, data = full_data[which(full_data$exp_arm=="Abstention" & full_data$task!=7),])
table2_abstention_NA_cluster_vcov <- cluster.vcov(table2_abstention_NA, cluster=full_data$ResponseId[which(full_data$exp_arm=="Abstention" & full_data$task!=7)])
coeftest(table2_abstention_NA, table2_abstention_NA_cluster_vcov)
table2_abstention_NA_ses <- coeftest(table2_abstention_NA, table2_abstention_NA_cluster_vcov)[,2]
table2_abstention_NA_pvals <- coeftest(table2_abstention_NA, table2_abstention_NA_cluster_vcov)[,4]

table2_abstention_NA_medium <- lm(choice_bin_NA ~ info_env*news, data = full_data[which(full_data$exp_arm=="Abstention" & full_data$task!=7 & full_data$info_env!="Low"),])
table2_abstention_NA_medium_cluster_vcov <- cluster.vcov(table2_abstention_NA_medium, cluster=full_data$ResponseId[which(full_data$exp_arm=="Abstention" & full_data$task!=7 & full_data$info_env!="Low")])
coeftest(table2_abstention_NA_medium, table2_abstention_NA_medium_cluster_vcov)
table2_abstention_NA_medium_ses <- coeftest(table2_abstention_NA_medium, table2_abstention_NA_medium_cluster_vcov)[,2]
table2_abstention_NA_medium_pvals <- coeftest(table2_abstention_NA_medium, table2_abstention_NA_medium_cluster_vcov)[,4]

print(
  texreg(l=list(table2_forced_NA, table2_forced_NA_medium,
              table2_abstention_NA, table2_abstention_NA_medium),
       custom.coef.map = list("(Intercept)" = "Intercept",
                              "info_envModerate" = "Medium information",
                              "info_envHigh" = "High information",
                              "newsRecently accused of cheating on spouse" = "Cheating on spouse",
                              "newsRecently accused of leaking confidential information" = "Leaking confidential information",
                              "newsRecently accused of sexual harassment" = "Sexual harassment",
                              "newsRecently celebrated wedding anniversary" = "Wedding anniversary",
                              "newsRecently honored for public service" = "Honored public service",
                              "info_envModerate:newsRecently accused of cheating on spouse" = "Medium information x cheating on spouse",
                              "info_envModerate:newsRecently accused of leaking confidential information" = "Medium information x leaking confidential information",
                              "info_envModerate:newsRecently accused of sexual harassment" = "Medium information x sexual harassment",
                              "info_envModerate:newsRecently celebrated wedding anniversary" = "Medium information x wedding anniversary",
                              "info_envModerate:newsRecently honored for public service" = "Medium information x honored public service",
                              "info_envHigh:newsRecently accused of cheating on spouse" = "High information x cheating on spouse",
                              "info_envHigh:newsRecently accused of leaking confidential information" = "High information x leaking confidential information",
                              "info_envHigh:newsRecently accused of sexual harassment" = "High information x sexual harassment",
                              "info_envHigh:newsRecently celebrated wedding anniversary" = "High information x wedding anniversary",
                              "info_envHigh:newsRecently honored for public service" = "High information x honored public service"),
       custom.model.names = c("Ref. Cat. = Low", "Ref. Cat. = Medium", 
                              "Ref. Cat. = Low", "Ref. Cat. = Medium"),
       override.se = list(table2_forced_NA_ses, table2_forced_NA_medium_ses,
                          table2_abstention_NA_ses, table2_abstention_NA_medium_ses),
       override.pvalues = list(table2_forced_NA_pvals, table2_forced_NA_medium_pvals,
                               table2_abstention_NA_pvals, table2_abstention_NA_medium_pvals),
       caption = "Funck and McCabe (2021) Replication Results (Table 2), Alternative Outcome Coding",
       caption.above = TRUE,
       stars = 0.05,
       include.rsquared = FALSE,
       include.adjrs = FALSE)
  , file = "../tables/table_SM6.tex"
)

#############
# FIGURE SM.3
#############

# estimating full model with interaction for experimental arm and outcome that
# codes abstentions as NAs, then retrieving corresponding standard errors and
# p-values

table2_arm_NA <- lm(choice_bin_NA ~ info_env*news*exp_arm, data = full_data[which(full_data$task!=7)])
table2_arm_NA_cluster_vcov <- cluster.vcov(table2_arm_NA, cluster=full_data$ResponseId[which(full_data$task!=7)])
coeftest(table2_arm_NA, table2_arm_NA_cluster_vcov)
table2_arm_NA_ses <- coeftest(table2_arm_NA, table2_arm_NA_cluster_vcov)[,2]
table2_arm_NA_pvals <- coeftest(table2_arm_NA, table2_arm_NA_cluster_vcov)[,4]

# calculating the differences in the AMCEs between the low and moderate and high
# information environments across the abstention and forced choice arms,
# separately, then plotting these differences

diffs_from_low_info_NA <- confint(glht(table2_arm_NA,
                                       linfct = c("`info_envModerate:newsRecently accused of sexual harassment` = 0",
                                                  "`info_envModerate:newsRecently accused of leaking confidential information` = 0",
                                                  "`info_envModerate:newsRecently accused of cheating on spouse` = 0",
                                                  "`info_envModerate:newsRecently celebrated wedding anniversary` = 0",
                                                  "`info_envModerate:newsRecently honored for public service` = 0",
                                                  "`info_envHigh:newsRecently accused of sexual harassment` = 0",
                                                  "`info_envHigh:newsRecently accused of leaking confidential information` = 0",
                                                  "`info_envHigh:newsRecently accused of cheating on spouse` = 0",
                                                  "`info_envHigh:newsRecently celebrated wedding anniversary` = 0",
                                                  "`info_envHigh:newsRecently honored for public service` = 0",
                                                  "`info_envModerate:newsRecently accused of sexual harassment` + `info_envModerate:newsRecently accused of sexual harassment:exp_armAbstention` = 0",
                                                  "`info_envModerate:newsRecently accused of leaking confidential information` + `info_envModerate:newsRecently accused of leaking confidential information:exp_armAbstention` = 0",
                                                  "`info_envModerate:newsRecently accused of cheating on spouse` + `info_envModerate:newsRecently accused of cheating on spouse:exp_armAbstention` = 0",
                                                  "`info_envModerate:newsRecently celebrated wedding anniversary` + `info_envModerate:newsRecently celebrated wedding anniversary:exp_armAbstention` = 0",
                                                  "`info_envModerate:newsRecently honored for public service` + `info_envModerate:newsRecently honored for public service:exp_armAbstention` = 0",
                                                  "`info_envHigh:newsRecently accused of sexual harassment` + `info_envHigh:newsRecently accused of sexual harassment:exp_armAbstention` = 0",
                                                  "`info_envHigh:newsRecently accused of leaking confidential information` + `info_envHigh:newsRecently accused of leaking confidential information:exp_armAbstention` = 0",
                                                  "`info_envHigh:newsRecently accused of cheating on spouse` + `info_envHigh:newsRecently accused of cheating on spouse:exp_armAbstention` = 0",
                                                  "`info_envHigh:newsRecently celebrated wedding anniversary` + `info_envHigh:newsRecently celebrated wedding anniversary:exp_armAbstention` = 0",
                                                  "`info_envHigh:newsRecently honored for public service` + `info_envHigh:newsRecently honored for public service:exp_armAbstention` = 0")),
                                  calpha = 2.81)$confint # 20 hypothesis tests, 0.05/20=0.0025, corresponding critical value 2.81

pdf(file = "../figures/figure_SM3.pdf", family = "Times", height = 11, width=9)
layout(matrix(c(1,2), nrow=2, byrow=TRUE), heights = c(0.90,0.10))
par(mar=c(6.1,12,1.75,1))
plot(x=c(diffs_from_low_info_NA[11,1],
         diffs_from_low_info_NA[1,1],
         diffs_from_low_info_NA[12,1],
         diffs_from_low_info_NA[2,1],
         diffs_from_low_info_NA[13,1],
         diffs_from_low_info_NA[3,1],
         diffs_from_low_info_NA[14,1],
         diffs_from_low_info_NA[4,1],
         diffs_from_low_info_NA[15,1],
         diffs_from_low_info_NA[5,1],
         diffs_from_low_info_NA[16,1],
         diffs_from_low_info_NA[6,1],
         diffs_from_low_info_NA[17,1],
         diffs_from_low_info_NA[7,1],
         diffs_from_low_info_NA[18,1],
         diffs_from_low_info_NA[8,1],
         diffs_from_low_info_NA[19,1],
         diffs_from_low_info_NA[9,1],
         diffs_from_low_info_NA[20,1],
         diffs_from_low_info_NA[10,1]), 
     y=c(1.3, 1.7,
         2.3, 2.7,
         3.3, 3.7,
         4.3, 4.7,
         5.3, 5.7,
         7.3, 7.7,
         8.3, 8.7,
         9.3, 9.7,
         10.3, 10.7,
         11.3, 11.7),
     xlim=c(-0.25,0.25),
     xaxt="n",
     ylim=c(1, 12),
     pch=c(rep(c(17,19))),
     col=c(rep(c("gray50","black"))),
     
     tck=-.02,
     cex.axis=0.9,
     cex=1.25,
     ylab="",
     yaxt="n",
     xlab="",
     axes = FALSE,
     panel.first = c(abline(v=0,lwd=2, col="gray70",lty=2)))
segments(col=c(rep(c("gray50","black"),18)),
         x0=c(diffs_from_low_info_NA[11,2],
              diffs_from_low_info_NA[1,2],
              diffs_from_low_info_NA[12,2],
              diffs_from_low_info_NA[2,2],
              diffs_from_low_info_NA[13,2],
              diffs_from_low_info_NA[3,2],
              diffs_from_low_info_NA[14,2],
              diffs_from_low_info_NA[4,2],
              diffs_from_low_info_NA[15,2],
              diffs_from_low_info_NA[5,2],
              diffs_from_low_info_NA[16,2],
              diffs_from_low_info_NA[6,2],
              diffs_from_low_info_NA[17,2],
              diffs_from_low_info_NA[7,2],
              diffs_from_low_info_NA[18,2],
              diffs_from_low_info_NA[8,2],
              diffs_from_low_info_NA[19,2],
              diffs_from_low_info_NA[9,2],
              diffs_from_low_info_NA[20,2],
              diffs_from_low_info_NA[10,2]),
         x1=c(diffs_from_low_info_NA[11,3],
              diffs_from_low_info_NA[1,3],
              diffs_from_low_info_NA[12,3],
              diffs_from_low_info_NA[2,3],
              diffs_from_low_info_NA[13,3],
              diffs_from_low_info_NA[3,3],
              diffs_from_low_info_NA[14,3],
              diffs_from_low_info_NA[4,3],
              diffs_from_low_info_NA[15,3],
              diffs_from_low_info_NA[5,3],
              diffs_from_low_info_NA[16,3],
              diffs_from_low_info_NA[6,3],
              diffs_from_low_info_NA[17,3],
              diffs_from_low_info_NA[7,3],
              diffs_from_low_info_NA[18,3],
              diffs_from_low_info_NA[8,3],
              diffs_from_low_info_NA[19,3],
              diffs_from_low_info_NA[9,3],
              diffs_from_low_info_NA[20,3],
              diffs_from_low_info_NA[10,3]),
         y0=c(1.3, 1.7,
              2.3, 2.7,
              3.3, 3.7,
              4.3, 4.7,
              5.3, 5.7,
              7.3, 7.7,
              8.3, 8.7,
              9.3, 9.7,
              10.3, 10.7,
              11.3, 11.7
         ))
axis(1,at=c(seq(-0.25,0.25,0.05)), 
     labels = c("-25%", NA, NA, NA, NA, "0%", NA, NA, NA, NA, "25%"), cex.axis = 1.5)
par(mgp=c(0,8,0))
axis(2,at=c(seq(1.5,5.5,1),
            seq(7.5,11.5,1)),
     las=2,
     labels=rep(c("Sexual Harassment",
                  "Leaking Confidential\nInformation",
                  "Cheating on Spouse",
                  "Celebrated Wedding\nAnniversary",
                  "Honored for Public\nService"),2),
     tck=0,
     lwd = 0,
     line = 1,
     cex.axis=1, hadj=0)
brackets(x1=-.26, x2=-.26, y1=1.1, y2=5.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Moderate Information")), side = 2, at=3.5, line = 10.2, cex = 1.75)
brackets(x1=-.26, x2=-.26, y1=7.1, y2=11.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("High Information")), side = 2, at=9.5, line = 10.2, cex = 1.75)
par(mgp=c(0,11,0))
mtext("Difference from\nLow Information AMCE", side = 1, at=0, line = 4, cex = 1.75)
par(mar=c(0,0,0,0))
plot(0,0, type="n", axes=FALSE, xlab="", ylab="")
legend(x=-0.06, y=0.8, legend=c("Forced Choice",
                                "Abstention Option"), lty=1, lwd=1, 
       pch = c(19,17),
       col = c("black", "grey50"), cex = 1.50)
dev.off()

#############
# FIGURE SM.4
#############

# to assess what task characteristcs might prompt abstention, we need a long
# data frame that includes information about both candidates in each task.  the
# following code is similar to the data cleaning and organizing at the top of this
# script, but instead of creating a separate data frame for each profile-task
# dyad, it creates a separate data frame for each task

long_data_task1 <- subset(data, select = c(ResponseId, party1_A, news1_A, party1_B, news1_B,  info_env1, choice1))
long_data_task1$task <- 1
colnames(long_data_task1) <- c("ResponseId", "partyA", "newsA", "partyB", "newsB",
                          "info_env", "choice", "task")

long_data_task2 <- subset(data, select = c(ResponseId, party2_A, news2_A, party2_B, news2_B,  info_env2, choice2))
long_data_task2$task <- 2
colnames(long_data_task2) <- c("ResponseId", "partyA", "newsA", "partyB", "newsB",
                          "info_env", "choice", "task")

long_data_task3 <- subset(data, select = c(ResponseId, party3_A, news3_A, party3_B, news3_B,  info_env3, choice3))
long_data_task3$task <- 3
colnames(long_data_task3) <- c("ResponseId", "partyA", "newsA", "partyB", "newsB",
                          "info_env", "choice", "task")

long_data_task4 <- subset(data, select = c(ResponseId, party4_A, news4_A, party4_B, news4_B,  info_env4, choice4))
long_data_task4$task <- 4
colnames(long_data_task4) <- c("ResponseId", "partyA", "newsA", "partyB", "newsB",
                          "info_env", "choice", "task")

long_data_task5 <- subset(data, select = c(ResponseId, party5_A, news5_A, party5_B, news5_B,  info_env5, choice5))
long_data_task5$task <- 5
colnames(long_data_task5) <- c("ResponseId", "partyA", "newsA", "partyB", "newsB",
                          "info_env", "choice", "task")

long_data_task6 <- subset(data, select = c(ResponseId, party6_A, news6_A, party6_B, news6_B,  info_env6, choice6))
long_data_task6$task <- 6
colnames(long_data_task6) <- c("ResponseId", "partyA", "newsA", "partyB", "newsB",
                          "info_env", "choice", "task")

long_data_task <- rbind(long_data_task1, long_data_task2, long_data_task3, long_data_task4, long_data_task5, 
                   long_data_task6)

full_data_task <- merge(demo_data, long_data_task, by = "ResponseId")

# discerning whether the respondent chose a candidate or abstained
full_data_task$choice_type <- ifelse(full_data_task$choice=="Candidate A" | full_data_task$choice=="Candidate B", "Candidate", NA)
full_data_task$choice_type <- ifelse(full_data_task$choice=="Neither candidate", "Abstention", full_data_task$choice_type)

# coding whether the respondent abstained
full_data_task$abstention <- ifelse(full_data_task$choice_type=="Abstention", 1, NA)
full_data_task$abstention <- ifelse(full_data_task$choice_type=="Candidate", 0, full_data_task$abstention)

# for each candidate, coding whether the candidate had a recent news item
# indicative of a scandal
full_data_task$candA_scandal <- ifelse(full_data_task$newsA=="Recently accused of cheating on spouse" |
                                      full_data_task$newsA=="Recently accused of sexual harassment" |
                                      full_data_task$newsA=="Recently accused of leaking confidential information", 
                                  1, NA)
full_data_task$candA_scandal <- ifelse(full_data_task$newsA=="Recently honored for public service" |
                                      full_data_task$newsA=="Recently celebrated wedding anniversary" |
                                      full_data_task$newsA=="No recent news", 
                                  0, full_data_task$candA_scandal)

full_data_task$candB_scandal <- ifelse(full_data_task$newsB=="Recently accused of cheating on spouse" |
                                      full_data_task$newsB=="Recently accused of sexual harassment" |
                                      full_data_task$newsB=="Recently accused of leaking confidential information", 
                                  1, NA)
full_data_task$candB_scandal <- ifelse(full_data_task$newsB=="Recently honored for public service" |
                                      full_data_task$newsB=="Recently celebrated wedding anniversary" |
                                      full_data_task$newsB=="No recent news", 
                                  0, full_data_task$candB_scandal)

# using respondents' PID to determine if each candidate in the task was a copartisan
full_data_task$copartA <- ifelse(full_data_task$partyA=="Democrat" & full_data_task$pid7>=5, 1, NA)
full_data_task$copartA <- ifelse(full_data_task$partyA=="Democrat" & full_data_task$pid7<=3, 0, full_data_task$copartA)
full_data_task$copartA <- ifelse(full_data_task$partyA=="Republican" & full_data_task$pid7<=3, 1, full_data_task$copartA)
full_data_task$copartA <- ifelse(full_data_task$partyA=="Republican" & full_data_task$pid7>=5, 0, full_data_task$copartA)

full_data_task$copartB <- ifelse(full_data_task$partyB=="Democrat" & full_data_task$pid7>=5, 1, NA)
full_data_task$copartB <- ifelse(full_data_task$partyB=="Democrat" & full_data_task$pid7<=3, 0, full_data_task$copartB)
full_data_task$copartB <- ifelse(full_data_task$partyB=="Republican" & full_data_task$pid7<=3, 1, full_data_task$copartB)
full_data_task$copartB <- ifelse(full_data_task$partyB=="Republican" & full_data_task$pid7>=5, 0, full_data_task$copartB)

# using information on partisanship and scandal to code whether each task had
# no scandalized candidates, only a copartisan scandalized candidate,
# only a non-copartisan scandalized candidate, or both candidates were scandalized
full_data_task$part_scandal <- ifelse(((full_data_task$copartA==1 & full_data_task$candA_scandal==1) &
                                      (full_data_task$copartB==0 & full_data_task$candB_scandal==1)) |
                                     ((full_data_task$copartA==0 & full_data_task$candA_scandal==1) &
                                          (full_data_task$copartB==1 & full_data_task$candB_scandal==1)), 
                                 "Both_scandal", NA)

full_data_task$part_scandal <- ifelse(((full_data_task$copartA==1 & full_data_task$candA_scandal==1) &
                                      (full_data_task$copartB==0 & full_data_task$candB_scandal==0)) |
                                     ((full_data_task$copartA==0 & full_data_task$candA_scandal==0) &
                                          (full_data_task$copartB==1 & full_data_task$candB_scandal==1)), 
                                 "Copart_scandal", full_data_task$part_scandal)

full_data_task$part_scandal <- ifelse(((full_data_task$copartA==1 & full_data_task$candA_scandal==0) &
                                      (full_data_task$copartB==0 & full_data_task$candB_scandal==1)) |
                                     ((full_data_task$copartA==0 & full_data_task$candA_scandal==1) &
                                          (full_data_task$copartB==1 & full_data_task$candB_scandal==0)), 
                                 "Noncopart_scandal", full_data_task$part_scandal)

full_data_task$part_scandal <- ifelse(((full_data_task$copartA==1 & full_data_task$candA_scandal==0) &
                                      (full_data_task$copartB==0 & full_data_task$candB_scandal==0)) |
                                     ((full_data_task$copartA==0 & full_data_task$candA_scandal==0) &
                                          (full_data_task$copartB==1 & full_data_task$candB_scandal==0)), 
                                 "No_scandal", full_data_task$part_scandal)

# specifying the ordering of the levels for information environment and scandal
full_data_task$info_env <- factor(full_data_task$info_env, levels = c("Low", "Moderate", "High"))
full_data_task$part_scandal <- factor(full_data_task$part_scandal, levels = c("No_scandal", "Copart_scandal", "Noncopart_scandal",
                                                                    "Both_scandal"))

# modeling abstention at the task level as a function of the prescence of scandal
# for partisan/non-copartisan candidates and the information environment
why_abstention <- lm(abstention ~ part_scandal*info_env, data = full_data_task[which(full_data_task$exp_arm=="Abstention"),])
summary(why_abstention)

# formal hypothesis tests for relationship between abstention and scandal
# both scandal low vs. both scandal high
summary(glht(why_abstention, linfct = "info_envHigh + part_scandalBoth_scandal:info_envHigh = 0"))
# both scandal low vs. both scandal medium
summary(glht(why_abstention, linfct = "info_envModerate + part_scandalBoth_scandal:info_envModerate = 0"))
# both scandal medium vs. both scandal high
summary(glht(why_abstention, linfct = "info_envModerate + part_scandalBoth_scandal:info_envModerate -
                                        info_envHigh - part_scandalBoth_scandal:info_envHigh= 0"))

# copart scandal low vs. copart scandal high
summary(glht(why_abstention, linfct = "info_envHigh + part_scandalCopart_scandal:info_envHigh = 0"))
# copart scandal low vs. copart scandal medium
summary(glht(why_abstention, linfct = "info_envModerate + part_scandalCopart_scandal:info_envModerate = 0"))
# copart scandal medium vs. copart scandal high
summary(glht(why_abstention, linfct = "info_envModerate + part_scandalCopart_scandal:info_envModerate -
                                        info_envHigh - part_scandalCopart_scandal:info_envHigh= 0"))

why_abstention_predictions <- predict(why_abstention, newdata = data.frame("part_scandal"=rep(c("No_scandal",
                                                                                                "Noncopart_scandal",
                                                                                                "Copart_scandal",
                                                                                                "Both_scandal"),3),
                                                                           "info_env"=c(rep("Low",4),
                                                                                        rep("Moderate",4),
                                                                                        rep("High",4))),
                                      se.fit = TRUE)

pdf(file = "../figures/figure_SM4.pdf", family = "Times", height = 8, width=8)
par(mar=c(6.1,14,2,1))
plot(x=c(why_abstention_predictions$fit[1:12]), 
     y=c(1:4, 6:9, 11:14),
     xlim=c(0,0.50),
     xaxt="n",
     ylim=c(1, 14),
     pch=19,
     col="black",
     tck=-.02,
     cex.axis=0.9,
     cex=1.25,
     ylab="",
     yaxt="n",
     xlab="",
     axes = FALSE)
segments(col="black",
         x0=c(why_abstention_predictions$fit[1:12]-1.96*why_abstention_predictions$se.fit[1:12]),
         x1=c(why_abstention_predictions$fit[1:12]+1.96*why_abstention_predictions$se.fit[1:12]),
         y0=c(1:4, 6:9, 11:14))
axis(1,at=c(0.00, 0.10, 0.20, 0.30, 0.40, 0.50), 
     labels = c("0.00", "0.10", "0.20", "0.30", "0.40", "0.50"), cex.axis = 1.5)
par(mgp=c(0,8,0))
axis(2,at=c(1:4, 6:9, 11:14),
     las=2,
     labels=rep(c("No_scandal",
                  "Non-copart scandal",
                  "Copart_scandal",
                  "Both_scandal"),3),
     tck=0,
     lwd = 0,
     line = 1,
     cex.axis=1, hadj=0)
brackets(x1=-.02, x2=-.02, y1=0.51, y2=4.49, type=4, ticks = NA, xpd=TRUE, h=0.01)
mtext(expression(~bold("Low Information")), side = 2, at=2.5, line = 10, cex = 1.25)
brackets(x1=-.02, x2=-.02, y1=5.51, y2=9.49, type=4, ticks = NA, xpd=TRUE, h=0.01)
mtext(expression(~bold("Moderate Information")), side = 2, at=7.5, line = 10, cex = 1.25)
brackets(x1=-.02, x2=-.02, y1=10.51, y2=14.49, type=4, ticks = NA, xpd=TRUE, h=0.01)
mtext(expression(~bold("High Information")), side = 2, at=12.5, line = 10, cex = 1.25)
par(mgp=c(0,11,0))
mtext("Probability of Abstention", side = 1, at=0.15, line = 4, cex = 1.75)
dev.off()

############
# TABLE SM.7
############

# recovering other demographic information from the original data frame
data$age <- ifelse(data$Q2=="18-34", 1, NA)
data$age <- ifelse(data$Q2=="35-44", 2, data$age)
data$age <- ifelse(data$Q2=="45-64", 3, data$age)
data$age <- ifelse(data$Q2=="65+", 4, data$age)

data$female <- ifelse(data$Q4=="Female", 1, NA)
data$female <- ifelse(data$Q4=="Male", 0, data$female)
data$female <- ifelse(data$Q4=="Other (Please specify)", 0, data$female)

data$education <- ifelse(data$Q6=="Some high school, or less", 1, NA)
data$education <- ifelse(data$Q6=="High school graduate or GED", 2, data$education)
data$education <- ifelse(data$Q6=="Some college, no 4-year degree", 3, data$education)
data$education <- ifelse(data$Q6=="College graduate", 4, data$education)
data$education <- ifelse(data$Q6=="Post-graduate degree", 5, data$education)

data$black <- ifelse(data$Q7=="Black or African-American", 1, NA)
data$black <- ifelse(data$Q7!="Black or African-American" & data$Q7!="", 0, data$black)

data$asian <- ifelse(data$Q7=="Asian", 1, NA)
data$asian <- ifelse(data$Q7!="Asian" & data$Q7!="", 0, data$asian)

data$otherrace <- ifelse(data$Q7=="American Indian or Alaska Native", 1, NA)
data$otherrace <- ifelse(data$Q7=="Native Hawaiian or other Pacific Islander", 1, data$otherrace)
data$otherrace <- ifelse(data$Q7=="Other (Please specify)", 1, data$otherrace)
data$otherrace <- ifelse(data$Q7!="Other (Please specify)" & 
                             data$Q7!="Native Hawaiian or other Pacific Islander" & 
                             data$Q7!="American Indian or Alaska Native" & 
                             data$Q7!="", 0, data$otherrace)

data$hispanic <- ifelse(data$Q8=="Yes", 1, NA)
data$hispanic <- ifelse(data$Q8=="No", 0, data$hispanic)

data$income <- ifelse(data$Q9=="Less than $20,000", 1, NA)
data$income <- ifelse(data$Q9=="$20,000-$29,999", 2, data$income)
data$income <- ifelse(data$Q9=="$30,000-$39,999", 3, data$income)
data$income <- ifelse(data$Q9=="$40,000-$49,999", 4, data$income)
data$income <- ifelse(data$Q9=="$50,000-$59,999", 5, data$income)
data$income <- ifelse(data$Q9=="$60,000-$69,999", 6, data$income)
data$income <- ifelse(data$Q9=="$70,000-$79,999", 7, data$income)
data$income <- ifelse(data$Q9=="$80,000-$89,999", 8, data$income)
data$income <- ifelse(data$Q9=="$90,000-$99,999", 9, data$income)
data$income <- ifelse(data$Q9=="$100,000-$124,999", 10, data$income)
data$income <- ifelse(data$Q9=="$125,000-$149,999", 11, data$income)
data$income <- ifelse(data$Q9=="$150,000 or more", 12, data$income)

data$ideology <- ifelse(data$Q11=="Very conservative", 1, NA)
data$ideology <- ifelse(data$Q11=="Conservative", 2, data$ideology)
data$ideology <- ifelse(data$Q11=="Slightly conservative", 3, data$ideology)
data$ideology <- ifelse(data$Q11=="Moderate", 4, data$ideology)
data$ideology <- ifelse(data$Q11=="Slightly liberal", 5, data$ideology)
data$ideology <- ifelse(data$Q11=="Liberal", 6, data$ideology)
data$ideology <- ifelse(data$Q11=="Very liberal", 7, data$ideology)

data$pid3 <- ifelse(data$Q14=="Strong Republican", "Republican", NA)
data$pid3 <- ifelse(data$Q14=="Not a very strong Republican", "Republican", data$pid3)
data$pid3 <- ifelse(data$Q15=="Republican", "Republican", data$pid3)
data$pid3 <- ifelse(data$Q15=="Neither", "Independent", data$pid3)
data$pid3 <- ifelse(data$Q15=="Democratic", "Democrat", data$pid3)
data$pid3 <- ifelse(data$Q13=="Not a very strong Democrat", "Democrat", data$pid3)
data$pid3 <- ifelse(data$Q13=="Strong Democrat", "Democrat", data$pid3)
data$pid3 <- ifelse(data$Q12=="Other", "Other", data$pid3)

# subsetting to the new set of demographics we want to add to our main data frame
demo_data <- subset(data, select = c(ResponseId, age, female, 
                                     education, black, asian, otherrace, 
                                     hispanic, income, ideology, pid3))

# merging demos with the task-level data frame used for Figure SM.4
full_data_task_demos <- merge(demo_data, full_data_task, by = "ResponseId")

# estimating abstention as a function of demographic characteristics
why_abstention <- lm(abstention ~ age + female + education + black + asian + otherrace +
                         hispanic + income + ideology, 
                     data = full_data_task_demos[which(full_data_task_demos$exp_arm=="Abstention"),])
summary(why_abstention)

print(
  texreg(l=list(why_abstention),
       custom.coef.map = list("(Intercept)" = "Intercept",
                              "age" = "Age",
                              "female" = "Female",
                              "education" = "Education",
                              "asian" = "Asian",
                              "black" = "Black",
                              "otherrace" = "Other Non-White Race",
                              "hispanic" = "Hispanic",
                              "income" = "Income",
                              "ideology" = "Ideology"
       ),
       caption = "Demographic correlates of abstention in F-M replication",
       caption.above = TRUE,
       stars = 0.05,
       include.rsquared = FALSE,
       include.adjrs = FALSE)
  , file = "../tables/table_SM7.tex"
)

#################
# IN-TEXT RESULTS
#################

# PAGE 10 AND SM.26--HOW COMMON ARE ABSTENTIONS?

# across all tasks 
table(full_data_task$abstention[which(full_data_task$exp_arm=="Abstention")])
prop.table(table(full_data_task$abstention[which(full_data_task$exp_arm=="Abstention")]))

# how many respondents abstained at least once?
ever_abstain <- aggregate(full_data_task$abstention[!is.na(full_data_task$abstention)], 
                          by = list(full_data_task$ResponseId[!is.na(full_data_task$abstention)]), FUN = "sum")
length(ever_abstain$Group.1[which(ever_abstain$x!=0)])
length(ever_abstain$Group.1[which(ever_abstain$x==0)])
length(ever_abstain$Group.1[which(ever_abstain$x!=0)])/dim(ever_abstain)[1]

# FOOTNOTE 6 AND SM.32--HOW OFTEN DO RESPONDENTS SWITCH FROM CHOOSING TO 
# ABSTAINING WHEN ALLOWED TO DO SO?

# coding attribute-levels and choices for seventh task
data$party7_A <- ifelse(str_detect(data$FL_180_DO, "D-R") | str_detect(data$FL_186_DO, "D-R"), "Democrat", NA)
data$party7_A <- ifelse(str_detect(data$FL_180_DO, "R-D") | str_detect(data$FL_186_DO, "R-D"), "Republican", data$party7_A)
data$news7_A <- data$`G-1-1-1`
data$party7_B <- ifelse(str_detect(data$FL_180_DO, "D-R") | str_detect(data$FL_186_DO, "D-R"), "Republican", NA)
data$party7_B <- ifelse(str_detect(data$FL_180_DO, "R-D") | str_detect(data$FL_186_DO, "R-D"), "Democrat", data$party7_B)
data$news7_B <- data$`G-1-2-1`
data$info_env7 <- ifelse(str_detect(data$FL_180_DO, "High") | str_detect(data$FL_186_DO, "High"), "High", NA)
data$info_env7 <- ifelse(str_detect(data$FL_180_DO, "Moderate") | str_detect(data$FL_186_DO, "Moderate"), "Moderate", data$info_env7)
data$info_env7 <- ifelse(str_detect(data$FL_180_DO, "Low") | str_detect(data$FL_186_DO, "Low"), "Low", data$info_env7)
data$choice7 <- ifelse(data$exp_arm=="Abstention", paste0(data$Q86, data$Q92, data$Q96, data$Q104, data$Q100, data$Q108),
                       paste0(data$Q235, data$Q239, data$Q243, data$Q247, data$Q251, data$Q255))

# how many in the forced choice arm chose to abstain when they were allowed to do so?
table(data$choice1[which(data$choice1!="" & data$choice7!="" & data$exp_arm=="Forced")], 
      data$choice7[which(data$choice1!="" & data$choice7!="" & data$exp_arm=="Forced")])

prop.table(table(data$choice1[which(data$choice1!="" & data$choice7!="" & data$exp_arm=="Forced")], 
                 data$choice7[which(data$choice1!="" & data$choice7!="" & data$exp_arm=="Forced")]))
